Hands-on 3 - Processing and Visualising Flow DataR

Modified

December 2, 2023

15.1 Overview

Spatial interaction represent the flow of people, material, or information between locations in geographical space. It encompasses everything from freight shipments, energy flows, and the global trade in rare antiquities, to flight schedules, rush hour woes, and pedestrian foot traffic.

Each spatial interaction, as an analogy for a set of movements, is composed of a discrete origin/destination pair. Each pair can be represented as a cell in a matrix where rows are related to the locations (centroids) of origin, while columns are related to locations (centroids) of destination. Such a matrix is commonly known as an origin/destination matrix, or a spatial interaction matrix.

In this hands-on exercise, you will learn how to build an OD matrix by using Passenger Volume by Origin Destination Bus Stops data set downloaded from LTA DataMall. By the end of this hands-on exercise, you will be able:

  • to import and extract OD data for a selected time interval,

  • to import and save geospatial data (i.e. bus stops and mpsz) into sf tibble data frame objects,

  • to populate planning subzone code into bus stops sf tibble data frame,

  • to construct desire lines geospatial data from the OD data, and

  • to visualise passenger volume by origin and destination bus stops by using the desire lines data.

15.2 Getting Started

For the purpose of this exercise, four r packages will be used. They are:

sf for importing, integrating, processing and transforming geospatial data. tidyverse for importing, integrating, wrangling and visualising data. tmap for creating thematic maps.

pacman::p_load(tmap, sf, DT, performance, ggpubr, stplanr, tidyverse)

15.3

Preparing the Flow Data 15.3.1 Importing the OD data Firstly, we will import the Passenger Volume by Origin Destination Bus Stops data set downloaded from LTA DataMall by using read_csv() of readr package.

odbus <- read_csv("data/aspatial/origin_destination_bus_202310.csv")

Let use display the odbus tibble data table by using the code chunk below.

glimpse(odbus)
Rows: 5,694,297
Columns: 7
$ YEAR_MONTH          <chr> "2023-10", "2023-10", "2023-10", "2023-10", "2023-…
$ DAY_TYPE            <chr> "WEEKENDS/HOLIDAY", "WEEKDAY", "WEEKENDS/HOLIDAY",…
$ TIME_PER_HOUR       <dbl> 16, 16, 14, 14, 17, 17, 17, 7, 14, 14, 10, 20, 20,…
$ PT_TYPE             <chr> "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "…
$ ORIGIN_PT_CODE      <chr> "04168", "04168", "80119", "80119", "44069", "2028…
$ DESTINATION_PT_CODE <chr> "10051", "10051", "90079", "90079", "17229", "2014…
$ TOTAL_TRIPS         <dbl> 3, 5, 3, 5, 4, 1, 24, 2, 1, 7, 3, 2, 5, 1, 1, 1, 1…
odbus$ORIGIN_PT_CODE <- as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- as.factor(odbus$DESTINATION_PT_CODE)

15.3.2 Extracting the study data

For the purpose of this exercise, we will extract commuting flows on weekday and between 6 and 9 o’clock.

odbus6_9 <- odbus %>% 
  filter(DAY_TYPE == "WEEKDAY") %>% 
  filter(TIME_PER_HOUR >= 6 & TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE, DESTINATION_PT_CODE) %>% 
  summarise(TRIPS = sum(TOTAL_TRIPS))

write_rds(odbus6_9,"data/rds/odbus6_9.rds")

odbus6_9 <- read_rds("data/rds/odbus6_9.rds")

Table below shows the content of odbus6_9

datatable(odbus6_9)

15.4 Working with Geospatial Data

For the purpose of this exercise, two geospatial data will be used. They are:

  • BusStop: This data provides the location of bus stop as at last quarter of 2022.

  • MPSZ-2019: This data provides the sub-zone boundary of URA Master Plan 2019. Both data sets are in ESRI shapefile format.

15.4.1 Importing geospatial data

Two geospatial data will be used in this exercise, they are:

busstop <- st_read(dsn =  "data/geospatial", layer = "BusStop") %>%
  st_transform(crs = 3414)
Reading layer `BusStop' from data source 
  `C:\weipengten\ISSS624\Hands-on_Ex\Hands-on_Ex3\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21
mpsz <- st_read(dsn =  "data/geospatial", layer = "MPSZ-2019") %>% 
  st_transform(crs = 3414)
Reading layer `MPSZ-2019' from data source 
  `C:\weipengten\ISSS624\Hands-on_Ex\Hands-on_Ex3\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

Note st_read() function of sf package is used to import the shapefile into R as sf data frame. st_transform() function of sf package is used to transform the projection to crs 3414.

15.5 Geospatial data wrangling

15.5.1 Combining Busstop and mpsz

Code chunk below populates the planning subzone code (i.e. SUBZONE_C) of mpsz sf data frame into busstop sf data frame.

busstop_mpsz <- st_intersection(busstop, mpsz) %>% 
  select(BUS_STOP_N, SUBZONE_C) %>% 
  st_drop_geometry()

write_rds(busstop_mpsz,  "data/rds/busstop_mpsz.rds")

Next, we are going to append the planning subzone code from busstop_mpsz data frame onto odbus6_9 data frame.

od_data <- left_join(odbus6_9 , busstop_mpsz, by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>% 
  rename(ORIGIN_BS = ORIGIN_PT_CODE, ORIGIN_SZ = SUBZONE_C, DESTIN_BS = DESTINATION_PT_CODE)

Before continue, it is a good practice for us to check for duplicating records.

duplicate <- od_data %>% 
  group_by_all() %>% 
  filter(n()>1) %>%
  ungroup()
duplicate
# A tibble: 1,022 × 4
   ORIGIN_BS DESTIN_BS TRIPS ORIGIN_SZ
   <chr>     <fct>     <dbl> <chr>    
 1 22501     22009         1 JWSZ09   
 2 22501     22009         1 JWSZ09   
 3 22501     22451       167 JWSZ09   
 4 22501     22451       167 JWSZ09   
 5 22501     22469        28 JWSZ09   
 6 22501     22469        28 JWSZ09   
 7 22501     22479        20 JWSZ09   
 8 22501     22479        20 JWSZ09   
 9 22501     22509         4 JWSZ09   
10 22501     22509         4 JWSZ09   
# ℹ 1,012 more rows

If duplicated records are found, the code chunk below will be used to retain the unique records.

od_data <- unique(od_data)

It will be a good practice to confirm if the duplicating records issue has been addressed fully.

Next, we will update od_data data frame cwith the planning subzone codes.

od_data <- left_join(od_data , busstop_mpsz, by = c("DESTIN_BS" = "BUS_STOP_N"))

duplicate <- od_data %>% 
  group_by_all() %>%
  filter(n()>1) %>% 
  ungroup()
duplicate
# A tibble: 1,324 × 5
   ORIGIN_BS DESTIN_BS TRIPS ORIGIN_SZ SUBZONE_C
   <chr>     <chr>     <dbl> <chr>     <chr>    
 1 01013     51071         2 RCSZ10    CCSZ01   
 2 01013     51071         2 RCSZ10    CCSZ01   
 3 01112     51071        66 RCSZ10    CCSZ01   
 4 01112     51071        66 RCSZ10    CCSZ01   
 5 01112     53041         4 RCSZ10    BSSZ01   
 6 01112     53041         4 RCSZ10    BSSZ01   
 7 01121     51071         8 RCSZ04    CCSZ01   
 8 01121     51071         8 RCSZ04    CCSZ01   
 9 01121     82221         1 RCSZ04    GLSZ05   
10 01121     82221         1 RCSZ04    GLSZ05   
# ℹ 1,314 more rows
od_data <- unique(od_data)

od_data <- od_data %>% 
rename(DESTIN_SZ = SUBZONE_C) %>% 
  drop_na() %>% 
  group_by(ORIGIN_SZ, DESTIN_SZ) %>% 
  summarise(MORNING_PEAK = sum(TRIPS))


write_rds(od_data,  "data/rds/od_data.rds")

od_data <- read_rds( "data/rds/od_data.rds")

15.6 Visualising Spatial Interaction

In this section, you will learn how to prepare a desire line by using stplanr package.

15.6.1 Removing intra-zonal flows

We will not plot the intra-zonal flows. The code chunk below will be used to remove intra-zonal flows.

od_data1 <- od_data[od_data$ORIGIN_SZ!=od_data$DESTIN_SZ,]

15.6.2 Creating desire lines

In this code chunk below, od2line() of stplanr package is used to create the desire lines.

flowLine <- od2line(flow = od_data1, zones = mpsz, zone_code = "SUBZONE_C")

15.6.3 Visualising the desire lines

To visualise the resulting desire lines, the code chunk below is used.

tm_shape(mpsz) + tm_polygons() + flowLine %>%
  tm_shape() + tm_lines(lwd = "MORNING_PEAK", style = "quantile", scale = c(0.1, 1, 3, 5, 7, 10), n = 6, alpha = 0.3)

When the flow data are very messy and highly skewed like the one shown above, it is wiser to focus on selected flows, for example flow greater than or equal to 5000 as shown below.

tm_shape(mpsz) + tm_polygons() + flowLine %>%
filter(MORNING_PEAK >= 5000) %>% 
  tm_shape() + 
  tm_lines(lwd = "MORNING_PEAK", 
           style = "quantile", 
           scale = c(0.1, 1, 3, 5, 7, 10), 
           n = 6, 
           alpha = 0.3) 

```